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1. INTRODUCTION 

There has been a great deal of research work in modeling the inelastic 
deformation behavior of materials at the elevated temperature environment. 

It appears that a unified approach in the form of viscoplastic relations has 
been most popular for prediction of material responses. In this context, a 
number of viscoplastic material models have been published in the literature 
[e.g. 1-6]. The unified approach differs from the conventional creep and 
plasticity theory in that both the creep and plastic deformations, or alter- 
nately termed inelastic deformations, are treated as time-dependent quanti- 
ties. Based on the experimental and theoretical studies performed by vari- 
ous investigators [3,4, 7-9], it is known that viscoplastic constitutive re- 
lations, in principle, are capable of predicting material responses at high 
temperatures such as cyclic plasticity, rate sensitivity, long-term creep 
deformations, strain-hardening or softening, etc. The degree of success of 
a constitutive relationship varies depending on the extent of parameters 
considered in or mathematical sophistication of a specific model. 

Although most of viscoplastic models give improved material response 
predictions over the classical approach, the associated constitutive differ- 
ential equations have stiff regimes which present numerical difficulties in 
time-dependent structural analysis. The numerical difficulty is indeed an 
important concern when the viscoplastic relations are applied to large scale 
finite element structural analysis. 
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In finite element analysis for viscoplastic materials, two issues of 
primary concern in connection with the associated material nonlinearity: 

1) solution convergence in solving the global (incremental) equilibrium 
equations, 2) integration of the constitutive rate equations at the local 
material points (or element integration points). Numerically, these two 
issues are inter-related. On the one hand, global equilibrium can not be 
achieved if the stresses calculated at local material points are grossly 
inaccurate. On the other hand, the constitutive relations and stresses are 
not representative to the material if the strains computed from the nodal 
displacements are in error. 

In view of the above discussion, we have therefore investigated a com- 
bined global/local incrementing scheme for the finite element analysis of 
viscoplastic materials. 

2. GLOBAL INCREMENTING 

Due to the material nonl inearity, a viscoplastic problem is effectively 
formulated by an incremental approach, in which the finite element equili- 
brium equations can be linearized. In order to solve these equations suc- 
cessfully, the analyst must be able to specify "appropriate" load steps. 

If the loading increments are too large, the solution may not converge, or 
it is far from being accurate. Alternatively, if the loading increments are 
very small, the computation cost will become prohibitively high. Therefore, 
it is desirable to implement an automatic incrementing procedure in which 
the selection of (global) load steps can be made by the program, rather than 
the analyst. 

Use of automatic load stepping for solving nonlinear problems is not 
new and most of the applications were concentrated at time independent 
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problems [10-11]. Herein, we adopted similar concepts for the solution of 
time dependent viscoplastic problems. The procedure involves two major 
steps: 1) initiation of incremental solution, and 2) selection of subsequent 
load increments. Each of the two steps is briefly outlined below. 

(1) Initiation of incremental solution - the solution begins with a speci- 
fied load vector, i.e. 


Rl = o • R (1) 

where a = a load factor, < 1. 

R = a reference load vector. 

With the above load vector, solution will proceed with equilibrium 
iterations. When the number of iterations reaches four and the solu- 
tion has not yet converged, an estimate is made to project the number 
of iterations required according to 

n = i + *n(DT0L/di)/(fcn dj - *n di_j) (2) 

where i = number of iterations already performed. 

DTOL = iteration tolerance for displacements, 
di = ratio between the incremental displacement norm of 
the i-th iteration and total displacement norm. 

- 1 4 «i i / 1 «i i 

If n is greater than a maximum number of iteration cycles allowed, then 
a new load vector is set to be RNEW = tRj, t<1. 
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(2) Subsequent load increments - The load increments, subsequent to the 
first step, are determined on the basis of a constant arc length 
method [12-14]. In this method, let the current load vector be 


R i+1 = M+l R (3) 

os 

where X-j+i = a load parameter corresponding to the (i+l)-th iteration 
= Xi + d X i+ i (4) 

and d x^+i is calculated from a quadratic algebraic equations [13]. 

3. LOCAL INCREMENTING 

Once the global load increment is determined from the method outlined 
in the above, a sub-incrementing method is incorporated at the material 
point level to integrate the rate constitutive equation. For the purpose of 
discussion, the viscoplastic constitutive equations are written in the form 

y * f (y» t) (5) 

/V os rv 

where y represents the vector of stress, inelastic strain and state vari- 
ables, and f is a vector of nonlinear functions. To integrate the preceding 
equations, we have developed an automatic procedure based on the variable- 
step Runge-Kutta (R-K) method. In this method, the global time increment At 
is divided into a number of sub-increments, i.e. h = At/n. Corresponding to 

h, the vector y for iteration (i+1) is evaluated by the 4th order and 5th 

(4) (5) 

order R-K formulas, respectively, i.e. y-j+i and y-j+i. Then an error can be 

os 

estimated from 

. | | (5) (4) 

E st - II yi+i - yi+i I l/ h (6) 

os OS 
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The calculated error must be within the following tolerance 

^st * 6 1 1 yi 1 1/ At (7) 

If the above condition is violated, a revised sub-increment h' is then 
obtained from 


h 1 s t h 


E s t * At 

The foregoing procedure is repeated until the criterion in Eq. (7) is 
satisfied. 


( 8 ) 


4. NUMERICAL EXAMP LE 

Several problems have been analyzed using the procedure outlined in 
the preceding sections. Presented herein is a thick walled cylinder sub- 
jected to an internal pressure, varying linearly from 0 to 14.6 psi for 
t e [0,40sec.] The cylinder material is assumed to be 2-1/4 Cr-Mo common 
steel at 811* k and Robinson's viscoplastic model is adopted. For finite 
element analysis, five 4-noded axisymmetric elements are used. 

The analysis was performed by using four different combinations of 
numerical algorithms: 

1) Automatic global and local incrementing (G + L) 

2) Automatic global incrementing (G) with constant local steps 
h = At/n s , n s = 2,4, and 8 

3) Automatic local sub- incrementing (L) with constant global steps, 

N • 5, 10, 16, and 20. 

4) Constant global and local steps. 

Summarized in Table 1 are the algorithm details, CPU time on IBM-3033 corn- 
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puter, and radial displacement at the outer surface of the cylinder. It is 
seen that with the automatic global/local incrementing algorithm, the lowest 
CPU time was consumed. Convergence difficulty was experienced if only the 
automatic global or local incrementing scheme was optioned unless the number 
of solution steps or the number of sub- increments is significantly in- 
creased. Shown in Figure 1 is the load vs. radial displacement at the outer 
surface of the cylinder calculated from two different algorithms, i.e. auto- 
matic global/local incrementing and constant global stepping with automatic 
local incrementing. Both algorithms gave almost identical results. 

5. CONCLUSION- 

Presented in this paper is a global/local time incrementing scheme for 
viscoplastic analysis of structures. The scheme is very efficient and use- 
ful for conducting large scale nonlinear finite element analysis involving 
viscoplastic materials. 
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Table 1 A COMPARISON OF DIFFERENT SOLUTION ALGORITHMS 
FOR A THICK-WALLED CYLINDER 


Case 

No. 

Option 

Global 

Steps 

■Hi 

mm 

1 

G + L 

16 

-- 

11 

0.1851 

2A 

G 

16 

8 

100 

0.1851 

2B 

G 

16 

4 

57 

0.1851 

2C 

G 

-- 

2 

solution diverged, (note 5) 

3A 

L 

20 

-- 

35 

0.1851 

3B 

L 

16 

-- 

36 

0.1856 

3C 

L 

10 

— 

45 

0.1879 

3D 

L 

5 

— 

solution diverged, (note 6) 

4A 

N 

20 

4 

37 

0.1851 

4B 

N 

20 

2 

solution diverged, (note 5) 


Note: 


1. G+L - both global and local automatic incrementing. 

2. L - local automatic incrementing only. 

3. G global automatic incrementing only. 

4. N - manual incrementing 

5. In cases 2C and 4B, solution diverged at steps 6 and 3, respectively, 
because the values of material state variables are out of bound. 

6. In case 3D, solution diverged at step 5 because out-of-balance load 
was greater than incremental load. 

7. U 0 is the radial displacement at outer surface of the cylinder. 
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